Tarea3

Minimos cuadrados

Plantear el problema de minimos cuadrados como un problema de minimizacion:

Tenemos que \(Y=X\beta + \epsilon\), por lo que donde \(\epsilon\) es el error, por lo que podemos plantear el problema como la minimizaci´ón de este error \[min_{ } ||\epsilon||^2 \] que es equivalente a \[min_{\beta} ||Y-X\beta||^2\] Para encontrar el estimador \(\hat{\beta}\) es necesario derivar e igualar a 0: Derivando respecto a \(\beta\): \[-2X^t (Y-X\hat{\beta})=0\] \[2X^tY=2X^tX\hat{\beta}\] \[\hat{\beta}=(X^tX)^{-1}X^tY\] Es el \(\hat{\beta}\) que soluciona el problema de minimización.

Con esto podemos ver que para cada: \[\hat{\beta}_{i}=\frac{\sum_{i}x_{i}y_{i}}{\sum_{i}x_{i}^2}\] La cuál es una composición lineal de los parametros. Desde el planteamiento estamos buscando que el resultado sea lineal en los parametros, no en los datos por lo que este metodo serviria para poder aproximar polinomios de la orma \(Y=X^2\)

¿En qué se relaciona esto con un problema de proyección?

Recordemos que al encontrar la proyección de un vector \(Proy_{X}(y)= \frac{<x,y>}{||x||^2}y\) el error de proyeccion \(Y- Proy_{x}(Y)\) será ortogonal a la proyección. En el caso del problema de minimizar el error, nosotros estamos proyectando la variable dependiente Y sobre el espacio formado por nuestros datos X. y buscamos las combinaciones sobre X que hagan este error mas pequeño.

Viendolo como el teorema de Pitagoras el cu´ál nos dice que el cuadrado de la hipotenisa de un triangulo es igual a la suma de los otros dos lados \(h^2=x^2+y^2\) si el angulo formado entre \(x\) y \(y\) es de 90 grados, en este caso \(H=Y\) las variables dependientes, \(x=Proy_{x}(Y)\) y \(y= error\), entonces la manera en la que se puede cumplir el teorema es que el error sea ortogonal a la proyeccion (angulo de 90 grados) sino el error seráa mas grande.

¿Que logramos al agregar una columna de unos en la matriz X?

Al definir una columna de unos estamos cambiando un poco el problema. Antes de agregar la columna: \[y_{j}=\beta_{1}x_{1j}+...+\beta_{n}x_{nj}\] Al definir la columa de unos estamos dando cierta holgura al modelo ya que al agregar un parametro \(\alpha\) no estamos forzando a que la aproximacion pase por el origen. \[y_{j}=\alpha+\beta_{1}x_{1j}+...+\beta_{n}x_{nj}\]

Plantear el problema de regresion como un problema de estadistica con errores….

¿Cual es la funcion de verosimilitud del problema anterior?

Planteando el problema como \[Y= X\beta +\epsilon\] con \[\epsilon \sim N(0,\sigma^2I_{n})\]

podemos probar que \[Y \sim N(X\beta,\sigma^2I_{n})\]

Demostracion: \[Y= X\beta +\epsilon\] \[\epsilon=Y-X\beta\]

1. Media: \[E(\epsilon)=E(Y)-E(X\beta)=0\] \[E(Y)=E(X\beta)\] \[E(Y)=X\beta\]
2. Varianza: \[Var(\epsilon)=Var(Y)+Var(X\beta)-2Cov(Y,X\beta)=\sigma^2I_{n}\] donde: \[Var(X\beta)=0\] y \[Cov(Y,X\beta)=E((Y-E(Y)(X\beta-E(X\beta)))\] \[=E((Y-X\beta)(X\beta-X\beta))=0\] entonces: \[Var(\epsilon)=Var(Y)=\sigma^2I_{n}\]

Ahora, falta demostrar que Y es Normal: Por propiedades de la Normal como \(\epsilon \sim N(0,\sigma^2I_{n})\) \[\epsilon + X\beta \sim N(X\beta, \sigma^2I_{n})\] (Esto era m´ás rapido pero bueno…)

Entonces con esto podemos escribir la funcion de verosimilitud de Y como: \[ L(\beta, \sigma^2)=\prod \frac{1}{\sqrt{2\pi\sigma^2}} exp(-\frac{(y-X\beta)^2}{2\sigma^2})\]

\[=(\frac {1}{\sqrt{2\pi\sigma^2}})^n exp (- \frac{1}{2\sigma} (y-X\beta)^t(y-X\beta))\]

A la cual le sacamos Logaritmo para despu´és poder buscar los parametros que maximicen la funcion: \[Log(L(\beta, \sigma^2))=-\frac{n}{2}ln(2\pi)-\frac{n}{2}ln(\sigma^2)-\frac{1}{2\sigma^2}(y-X\beta)^t(y-X\beta)\] Ahora hay que maximizar esa funcion sobre \(\beta\) y \(\sigma^2\) Respecto de \(\beta\) \[\frac{dL}{d\beta}=\frac{1}{\sigma^2}(y-X\beta)^tX=0\] \[X^t y=X^tX\beta\] \[\hat{\beta}=(X^tX)^{-1}X^ty\]

Respecto de \(\sigma^2\) \[\frac{dL}{d\sigma^2}= - \frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}(y-X\beta)^t(y-X\beta)=0\] \[\sigma^2=\frac{1}{n} (y-X\beta)^t(y-X\beta)\] que puede resolverse usando \(\hat{\beta}\) que encontramos en la solucion anterior.

Como pudimos ver el resultado de maxima verosimilitud \[\hat{\beta}=(X^tX)^{-1}X^ty\] Es igual al resultado de \(\hat{\beta}\) de minimos cuadrados

El teorema de Gauss Markov

El teorema de Gauss Markov plantea que un modelo de regresion lineal en el que se cumple que:

1. \(E(\epsilon)=0\)
2. Los errores no estan correlacionados entre si
3. \(Var(\epsilon_{i}=\sigma^2\) para todo i, los errores tienen la misma varianza

Entonces el mejor estimador lineal insengado de los oeficientes \(\beta\) es el resultante del problema de minimizar los errores al cuadrado. si existe.

Parte Aplicada

library(ggplot2)
library(plotly)
data("diamonds")
head(diamonds)

#Entonces, para hacer la regresion lineal tengo que agarrar los #valores numericos que son
diamonds_num<-diamonds[,c(1,5:10)]
reg_precio=lm(price~carat+depth+table+x+y+z, data=diamonds_num)
summary(reg_precio)

Call:
lm(formula = price ~ carat + depth + table + x + y + z, data = diamonds_num)

Residuals:
     Min       1Q   Median       3Q      Max 
-23878.2   -615.0    -50.7    347.9  12759.2 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 20849.316    447.562  46.584  < 2e-16 ***
carat       10686.309     63.201 169.085  < 2e-16 ***
depth        -203.154      5.504 -36.910  < 2e-16 ***
table        -102.446      3.084 -33.216  < 2e-16 ***
x           -1315.668     43.070 -30.547  < 2e-16 ***
y              66.322     25.523   2.599  0.00937 ** 
z              41.628     44.305   0.940  0.34744    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1497 on 53933 degrees of freedom
Multiple R-squared:  0.8592,    Adjusted R-squared:  0.8592 
F-statistic: 5.486e+04 on 6 and 53933 DF,  p-value: < 2.2e-16

Que tan bien ajusta el modelo depende de la correlacion que tenga cada una de las variables explicativas sobre el precio, esto lo podemos ver a trav´és de:

diamonds_num<- diamonds_num[,c(4,1,2,3,5,6,7)]
cor(diamonds_num)

           price      carat       depth      table           x
price  1.0000000 0.92159130 -0.01064740  0.1271339  0.88443516
carat  0.9215913 1.00000000  0.02822431  0.1816175  0.97509423
depth -0.0106474 0.02822431  1.00000000 -0.2957785 -0.02528925
table  0.1271339 0.18161755 -0.29577852  1.0000000  0.19534428
x      0.8844352 0.97509423 -0.02528925  0.1953443  1.00000000
y      0.8654209 0.95172220 -0.02934067  0.1837601  0.97470148
z      0.8612494 0.95338738  0.09492388  0.1509287  0.97077180
                y          z
price  0.86542090 0.86124944
carat  0.95172220 0.95338738
depth -0.02934067 0.09492388
table  0.18376015 0.15092869
x      0.97470148 0.97077180
y      1.00000000 0.95200572
z      0.95200572 1.00000000

Como podemos ver en la tabla de correlaciones, precio est´á correlacionada altamente con “carat”, “x”,“y”,“z” y tienen cierta correlacion con death y table, esto me dice que las variables escogidas sirven para explicar el precio. Podemos ver esta relacion de manera grafica a continuacion en el siguiente Diagrama de dispersion :

library(GGally)

package ‘GGally’ was built under R version 3.2.5replacing previous import by ‘utils::capture.output’ when loading ‘GGally’replacing previous import by ‘utils::head’ when loading ‘GGally’replacing previous import by ‘utils::installed.packages’ when loading ‘GGally’replacing previous import by ‘utils::str’ when loading ‘GGally’

library(ggplot2)
my_fn <- function(data, mapping, ...){
  p <- ggplot(data = data, mapping = mapping) + 
    geom_point() + 
    geom_smooth(method=loess, fill="red", color="red", ...) +
    geom_smooth(method=lm, fill="blue", color="blue", ...)
  p
}
g = ggpairs(diamonds_num,columns = 1:7, lower = list(continuous = my_fn))
g

 plot: [1,1] [=--------------------------------]  2% est: 0s 
 plot: [1,2] [=--------------------------------]  4% est: 5s 
 plot: [1,3] [==-------------------------------]  6% est: 5s 
 plot: [1,4] [===------------------------------]  8% est: 5s 
 plot: [1,5] [===------------------------------] 10% est: 5s 
 plot: [1,6] [====-----------------------------] 12% est: 5s 
 plot: [1,7] [=====----------------------------] 14% est: 4s 
 plot: [2,1] [=====----------------------------] 16% est: 4s 
 plot: [2,2] [======---------------------------] 18% est: 2m 
 plot: [2,3] [=======--------------------------] 20% est: 2m 
 plot: [2,4] [=======--------------------------] 22% est: 2m 
 plot: [2,5] [========-------------------------] 24% est: 1m 
 plot: [2,6] [=========------------------------] 27% est: 1m 
 plot: [2,7] [=========------------------------] 29% est: 1m 
 plot: [3,1] [==========-----------------------] 31% est: 1m 
 plot: [3,2] [===========----------------------] 33% est: 2m 
 plot: [3,3] [===========----------------------] 35% est: 3m 
 plot: [3,4] [============---------------------] 37% est: 2m 
 plot: [3,5] [=============--------------------] 39% est: 2m 
 plot: [3,6] [=============--------------------] 41% est: 2m 
 plot: [3,7] [==============-------------------] 43% est: 2m 
 plot: [4,1] [===============------------------] 45% est: 2m 
 plot: [4,2] [===============------------------] 47% est: 2m 
 plot: [4,3] [================-----------------] 49% est: 2m 
 plot: [4,4] [=================----------------] 51% est: 3m 
 plot: [4,5] [==================---------------] 53% est: 2m 
 plot: [4,6] [==================---------------] 55% est: 2m 
 plot: [4,7] [===================--------------] 57% est: 2m 
 plot: [5,1] [====================-------------] 59% est: 2m 
 plot: [5,2] [====================-------------] 61% est: 2m 
 plot: [5,3] [=====================------------] 63% est: 2m 
 plot: [5,4] [======================-----------] 65% est: 2m 
 plot: [5,5] [======================-----------] 67% est: 2m 
 plot: [5,6] [=======================----------] 69% est: 2m 
 plot: [5,7] [========================---------] 71% est: 2m 
 plot: [6,1] [========================---------] 73% est: 2m 
 plot: [6,2] [=========================--------] 76% est: 2m 
 plot: [6,3] [==========================-------] 78% est: 2m 
 plot: [6,4] [==========================-------] 80% est: 1m 
 plot: [6,5] [===========================------] 82% est: 1m 
 plot: [6,6] [============================-----] 84% est: 1m 
 plot: [6,7] [============================-----] 86% est: 1m 
 plot: [7,1] [=============================----] 88% est: 1m 
 plot: [7,2] [==============================---] 90% est:48s 
 plot: [7,3] [==============================---] 92% est:40s 
 plot: [7,4] [===============================--] 94% est:31s 
 plot: [7,5] [================================-] 96% est:22s 
 plot: [7,6] [================================-] 98% est:11s 
 plot: [7,7] [=================================]100% est: 0s 
                                                             

Graficamente así se ven los valores del modelo contra los observados. Podemos ver que no forman una linea de 45 grados por lo que posiblemente hay alg´ún problema dada la correlacion entre las variables explicativas. Esto tambien lo podemos ver en la tabla de correlaciones.

#Grafica de los valores del modelo contra los reales
library(plotly)
data<-as.data.frame(cbind(reg_precio$fitted.values,diamonds_num$price))
p <- plot_ly(data = data, x = ~data[,1], y = ~data[,2]) %>%
  layout(title="Valores ajustadosvs Valores observados")

replacing previous import by ‘shiny::includeHTML’ when loading ‘crosstalk’replacing previous import by ‘shiny::knit_print.shiny.tag’ when loading ‘crosstalk’replacing previous import by ‘shiny::code’ when loading ‘crosstalk’replacing previous import by ‘shiny::includeScript’ when loading ‘crosstalk’replacing previous import by ‘shiny::includeMarkdown’ when loading ‘crosstalk’replacing previous import by ‘shiny::tags’ when loading ‘crosstalk’replacing previous import by ‘shiny::is.singleton’ when loading ‘crosstalk’replacing previous import by ‘shiny::withTags’ when loading ‘crosstalk’replacing previous import by ‘shiny::img’ when loading ‘crosstalk’replacing previous import by ‘shiny::tagAppendAttributes’ when loading ‘crosstalk’replacing previous import by ‘shiny::knit_print.shiny.tag.list’ when loading ‘crosstalk’replacing previous import by ‘shiny::knit_print.html’ when loading ‘crosstalk’replacing previous import by ‘shiny::tagAppendChild’ when loading ‘crosstalk’replacing previous import by ‘shiny::includeCSS’ when loading ‘crosstalk’replacing previous import by ‘shiny::br’ when loading ‘crosstalk’replacing previous import by ‘shiny::singleton’ when loading ‘crosstalk’replacing previous import by ‘shiny::span’ when loading ‘crosstalk’replacing previous import by ‘shiny::a’ when loading ‘crosstalk’replacing previous import by ‘shiny::tagList’ when loading ‘crosstalk’replacing previous import by ‘shiny::strong’ when loading ‘crosstalk’replacing previous import by ‘shiny::tag’ when loading ‘crosstalk’replacing previous import by ‘shiny::p’ when loading ‘crosstalk’replacing previous import by ‘shiny::validateCssUnit’ when loading ‘crosstalk’replacing previous import by ‘shiny::HTML’ when loading ‘crosstalk’replacing previous import by ‘shiny::h1’ when loading ‘crosstalk’replacing previous import by ‘shiny::h2’ when loading ‘crosstalk’replacing previous import by ‘shiny::h3’ when loading ‘crosstalk’replacing previous import by ‘shiny::h4’ when loading ‘crosstalk’replacing previous import by ‘shiny::h5’ when loading ‘crosstalk’replacing previous import by ‘shiny::h6’ when loading ‘crosstalk’replacing previous import by ‘shiny::tagAppendChildren’ when loading ‘crosstalk’replacing previous import by ‘shiny::em’ when loading ‘crosstalk’replacing previous import by ‘shiny::div’ when loading ‘crosstalk’replacing previous import by ‘shiny::pre’ when loading ‘crosstalk’replacing previous import by ‘shiny::htmlTemplate’ when loading ‘crosstalk’replacing previous import by ‘shiny::suppressDependencies’ when loading ‘crosstalk’replacing previous import by ‘shiny::tagSetChildren’ when loading ‘crosstalk’replacing previous import by ‘shiny::includeText’ when loading ‘crosstalk’replacing previous import by ‘shiny::hr’ when loading ‘crosstalk’

p

No trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
No trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

Podemos ver que como el modelo tiene una \(R^2 ~ .8592\) por lo que podemos decir que el modelo explica 85% de la varianza con las variables que estamos usando pero hay algunos problemas entre ellas por lo que los estimadores deben de estar sesgados. nuestro error estandar residual \(\sigma=1497\) por lo que \(sigma^2=1497^2\)

¿Cual es el angulo entre Y y \(\hat{Y}\)

\[arcos(R^2)=\theta\] En este caso en particular \(R^2=.8592\) por lo que

acos(sqrt(.8592))

[1] 0.3846484

Haciendo una funcion de logverosimilitud

En caso de que podamos aproximar con una normal tenemos lo siguiente:

lik1<-function(par,X,y){
  beta<-par[1:ncol(X)]
  sigma2<-par[ncol(X)+1]
  n<-nrow(y)
  X<-as.matrix(X)
  beta<-as.matrix(beta)
  
  mu<-X%*%beta
  
  logl= -.5*n*log(2*pi) -.5*n*log(sigma2) - (sum((y-mu)**2)* .5/sigma2))

Error: unexpected ')' in:
"
  logl= -.5*n*log(2*pi) -.5*n*log(sigma2) - (sum((y-mu)**2)* .5/sigma2))"

En nuestro caso en particular:

optim(par=betasigma,lik1, X=X_reg, y=y_reg,method="BFGS")

NaNs producedNaNs producedNaNs produced

$par
[1]  27260.5528  10906.0130   -271.6052   -134.3417  -1586.6727    271.2380 140486.4120

$value
[1] 802048.1

$counts
function gradient 
     159      100 

$convergence
[1] 1

$message
NULL

reg_precio$coefficients

(Intercept)       carat       depth       table           x           y           z 
 20849.3164  10686.3091   -203.1541   -102.4457  -1315.6678     66.3216     41.6277 

Ahora, como z es “no significativa” vamos a sacarla del modelo y a correr todo con el nuevo conjunto de variables

summary(reg_red)

Call:
lm(formula = price ~ carat + depth + table + x + y, data = diamonds_red)

Residuals:
     Min       1Q   Median       3Q      Max 
-23872.1   -614.8    -50.5    347.5  12759.4 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 20702.947    419.575  49.343  < 2e-16 ***
carat       10686.707     63.199 169.095  < 2e-16 ***
depth        -200.718      4.855 -41.344  < 2e-16 ***
table        -102.490      3.084 -33.234  < 2e-16 ***
x           -1293.542     36.063 -35.869  < 2e-16 ***
y              69.575     25.287   2.751  0.00594 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1497 on 53934 degrees of freedom
Multiple R-squared:  0.8592,    Adjusted R-squared:  0.8592 
F-statistic: 6.583e+04 on 5 and 53934 DF,  p-value: < 2.2e-16

Ahora con esto vuelvo a definir mis variables a usar en la funcion de max verosimilitud

optim(par=betasigma,lik1, X=X_reg, y=y_reg, method = "BFGS")

NaNs producedNaNs producedNaNs produced

$par
[1]  27260.5528  10906.0130   -271.6052   -134.3417  -1586.6727    271.2380 140486.4120

$value
[1] 802048.1

$counts
function gradient 
     159      100 

$convergence
[1] 1

$message
NULL

Podemos ver como los estimadores obtenidos por maxima verosimilitud se parecen mas a los que tienen mayor significancia en el modelo de regresion lineal.